Astronomy & Astrophysics manuscript no. heflash'poplapoplll 


©ESO 2010 


March 19, 2010 





The core helium flash revisited 

III. From Pop I to Pop III stars 

M. Mocak 1 , S. W. Campbell 23 , E. Miiller 4 and K. Kifonidis 4 

1 Institut d' Astronomie et d' Astrophysique, Universite Libre de Bruxelles, CP 226, 1050 Brussels, Belgium 
e-mail: mmocakOulb . ac . be 

2 Departament de Ffsica i Enginyeria Nuclear, EUETIB, Universitat Politecnica de Catalunya, C./Comte d'Urgell 187, E-08036 
Barcelona, Spain 

e-mail: simon . w . campbell@upc . edu 

3 Centre for Stellar and Planetary Astrophysics, School of Mathematical Sciences, Monash University, Melbourne 3800, Australia 

4 Max-Planck-Institut fiir Astrophysik, Postfach 1312, 85741 Garching, Germany 



Received 

ABSTRACT 

Context. Degenerate ignition of helium in low-mass stars at the end of the red giant branch phase leads to dynamic convection in their 
helium cores. One-dimensional (ID) stellar modeling of this intrinsically multi-dimensional dynamic event is likely to be inadequate. 
Previous hydrodynamic simulations imply that the single convection zone in the helium core of metal-rich Pop I stars grows during 
the flash on a dynamic timescale. This may lead to hydrogen injection into the core, and a double convection zone structure as known 
from one-dimensional core helium flash simulations of low-mass Pop III stars. 

Aims. We perform hydrodynamic simulations of the core helium flash in two and three dimensions to better constrain the nature 
of these events. To this end we study the hydrodynamics of convection within the helium cores of a 1.25 M G metal-rich Pop I star 
(Z=0.02), and a 0.85 M Q metal-free Pop III star (Z=0) near the peak of the flash. These models possess single and double convection 
zones, respectively. 

Methods. We use ID stellar models of the core helium flash computed with state-of-the-art stellar evolution codes as initial models for 
our multidimensional hydrodynamic study, and simulate the evolution of these models with the Riemann solver based hydrodynamics 
code Herakles which integrates the Euler equations coupled with source terms corresponding to gravity and nuclear burning. 
Results. The hydrodynamic simulation of the Pop I model involving a single convection zone covers 27 hours of stellar evolution, 
while the first hydrodynamic simulations of a double convection zone, in the Pop III model, span 1.8 hours of stellar life. We find 
differences between the predictions of mixing length theory and our hydrodynamic simulations. The simulation of the single convec- 
tion zone in the Pop I model shows a strong growth of the size of the convection zone due to turbulent entrainment. Hence we predict 
that for the Pop I model a hydrogen injection phase (i.e., hydrogen injection into the helium core) will commence after about 23 days, 
which should eventually lead to a double convection zone structure known from ID stellar modeling of low-mass Pop III stars. Our 
two and three-dimensional hydrodynamic simulations of the double (Pop III) convection zone model show that the velocity field in 
the convection zones is different from that predicted by stellar evolutionary calculations. The simulations suggest that the double 
convection zone decays quickly, the flow eventually being dominated by internal gravity waves. 

Key words. Stars: evolution - hydrodynamics - convection - turbulent entrainment 



1. Introduction 

Runaway nuclear burning of helium in the core of low mass red 
giant stars leads to convective mixing and burning on dynamic 
time scales. One dimensional evolutionary simulations (which 
assume time scales much longer than the dynamical ones) may 
miss key features of this rapid phase that could have significant 
effects on the further evolution of the stars. Furthermore, ID hy- 
drodynamical simulations of this intrinsically multi-dimensional 
event is likely to be inadequate. 

Our previous hydrodynamic simulations (Mocak et al. 2008, 
2009) imply that a 1.25 M solar-like star can experience injec- 
tion of hydrogen into its helium core during the canonical core 
helium flash near its peak. Hydrogen injection results from the 
growth of the convection zone (which is sustained by helium 
burning) due to turbulent entrainment on a dynamic timescale 
(Meakin & Arnett 2007), and probably occurs for all low-mass 
Pop I stars, as the properties of their cores are similar at the peak 



of the core helium flash (Sweigart & Gross 1978). An obvious 
consequence of this scenario is that the convection zones are en- 
larged in these stars. Whether they fail to dredge up nuclear ash 
to the atmosphere shortly after the flash is still unclear. However, 
such a dredge up could explain the Al-Mg anticorrelation found 
in red giants at the tip of the red giant branch (Shetrone 1996a,b; 
Yong et al. 2006). In ID simulations one has to manipulate the 
properties of the core helium flash to achieve such a dredge up, 
e.g., by changing the ignition position of the helium in the core 
(Paczynski & Tremaine 1977) or by forcing inward mixing of 
hydrogen (Fujimoto et al. 1999). 



Canonical (ID) stellar evolution calculations predict hy- 
drogen injection during the core helium flash and subsequent 
dredge-up of nuclear ashes to the atmosphere only for Pop III 
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1 and extremely metal-poor (EMP; with intrinsic metallicities 
[Fe/H] < -4) stars. This is a promising scenario for explain- 
ing the peculiar abundances of carbon and nitrogen observed in 
Galactic EMP Halo stars (Campbell & Lattanzio 2008). If these 
stars are assumed to be polluted by accretion of CNO-rich inter- 
stellar matter they will possibly experience hydrogen injection 
but no subsequent dredge-up, because a large CNO metallicitiy 
(as compared to the intrinsic [Fe/H] metallicity) in the stellar en- 
velope influences the ignition site of the first major core helium 
flash, and thereby the occurrence of the dredge-up (Hollowell 
et al. 1990). The helium abundance adopted in the stellar models 
also seems to influence the process of hydrogen injection itself as 
shown by Schlattl et al. (2001), while the same authors find that 
the injection process seems to be independent of the assumed 
convection model. 

Stellar models with a higher intrinsic metallicity, i.e., [Fe/H] 
> -4, do not inject hydrogen into the helium core, and conse- 
quently there is also no dredge-up of CNO-rich nuclear ashes 
to the atmosphere (Fujimoto et al. 1990; Hollowell et al. 1990; 
Campbell & Lattanzio 2008). Whether this is the final answer re- 
mains unclear, however, as Fujimoto et al. (1999) with his semi- 
analytic study and a postulated hydrogen injection followed by 
a dredge-up could show that such a scenario can explain some 
peculiarities observed in the spectra of red-giant stars (related to 
CNO elements and 24 Mg) with metallicities as large as [Fe/H] 
<-l. 

There exist two main reasons why hydrogen injection 
episodes occur only in Pop III and EMP stars: (z) these stars 
possess a flatter entropy gradient in the hydrogen burning shell, 
and (ii) they ignite helium far off center, relatively close to the 
hydrogen-rich envelope (Fujimoto et al. 1990). However, Pop II 
and Pop I stars could also mix hydrogen into the helium core 
during the core helium flash 2 
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Fig. 1. Upper panel: Kippenhahn diagram of a stellar evolution- 
ary calculation during the core helium flash of a 0.85 M Pop III 
star with convection zones sustained by helium (He-rich) and 
hydrogen (H-rich) or CNO burning (grey shaded regions, except 
for the convective envelope). The border between the helium and 
hydrogen rich layers is given by the solid blue curve. The loca- 
tion of the initial model SC (model number 9120) - black vertical 
line. Lower panel: the temporal evolution of the helium (dotted- 
blue) and hydrogen (solid-red) luminosity as a function of time. 



- if the flash was more violent, and thus the helium convection 
zone wider (Despain & Scalo 1976; Despain 1981). This sce- 
nario is disfavored by the facts that the flash is less violent in 
stars with higher metallicity as less energy is needed to lift 
the degeneracy of the less massive cores (Sweigart & Gross 
1978), and that helium ignition occurs at smaller densities 
(Fujimoto et al. 1990), 

- or if the entropy gradient between the hydrogen and helium 
burning shell was sufficiently shallow (Iben 1976; Fujimoto 
1977). A small entropy gradient would allow the convective 
shell in the helium core to reach the hydrogen layers even 
though the flash itself would not be very violent. This sce- 
nario is also disfavored as solutions to the stellar structure 
equations are robust with many different groups getting very 
similar results i.e., no hydrogen injection (Fujimoto et al. 
1990; Hollowell et al. 1990; Campbell & Lattanzio 2008) 

- or if a growth of the helium convection zone through turbu- 
lent entrainment at the convective boundaries (Mocak et al. 
2008, 2009) could be sustained for a sufficiently long period 
of time. 



1 They are supposed to be the first stars in the Universe and seem to 
be extremely rare, as the most metal-poor star discovered up to now has 
a metallicity of [Fe/H] ~ -5.5 (Frebel et al. 2005). 

2 According to Schlattl et al. (2001) the occurrence of a hydrogen 
flash is favored by a higher electron degeneracy in the helium core, 
which leads to helium ignition closer to the hydrogen shell. 



A hydrogen injection phase also occurs in low-mass metal 
deficient stars on the asymptotic giant branch (AGB) at the be- 
ginning of the thermally pulsing stage (TPAGB) 3 . 

Hydrogen injection is found to occur in more massive stars 
(M > 1.3 Mo) with low metallicity during the TPAGB (Chieffi 
et al. 2001; Siess et al. 2002; Iwamoto et al. 2004), in "Late Hot 
Flasher" stars experiencing strong mass loss on the RGB (Brown 
et al. 2001; Cassisi et al. 2003), and in H-deficient post AGB 
(PAGB) stars. These events are referred to with various names 
in the literature. Here we use the nomenclature "dual flashes" 
(Campbell & Lattanzio 2008), since they all have in common 
simultaneous hydrogen and helium flashes. 

Dual flash events often lead to a splitting of the single he- 
lium convection zone (HeCZ) into two parts (double convection 
zone): one sustained by helium burning, and a second one by 
hydrogen burning via CNO cycles (Fig. 1). Double convection 
zones are structures which are commonly encountered in stellar 
models, but their hydrodynamic properties have so far only been 
studied for the oxygen and carbon burning shell of a 23 M Q star 
by Meakin & Arnett (2006). 

In the following we describe two-dimensional (2D) and 
three-dimensional (3D) hydrodynamic simulations of a helium 
core during the core helium flash with a single convection zone 
(Pop I; in 3D only) and a double convection zone (Pop III, in 

3 If hydrogen injection occurs at the tip of the red giant branch branch 
(RGB), it does not occur on the AGB - the star evolves like a normal 
thermally pulsating Pop I or II star (Schlattl et al. 2001) 
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2D and 3D), respectively. Previous studies have indicated that 
there is a strong interaction between the adjacent shells of a dou- 
ble convection zone by internal gravity waves (Meakin & Arnett 
2006). 

We introduce the stellar models used as input for our hy- 
drodynamic simulations in Sect. 2, briefly discuss the physics 
included in our simulations in Sect. 3, and give a short descrip- 
tion of our hydrodynamics code and the computational setup in 
Sect. 4. Subsequently, we present and compare the results of our 
2D and 3D hydrodynamic simulations in Sect. 5. In particular, 
we discuss turbulent entrainment at the convective boundaries 
for our single convection zone model, the temporal evolution of 
its kinetic energy density, and how our results compare with the 
predictions of mixing-length theory (MLT). We proceed simi- 
larly for our hydrodynamic double convection zone model, ex- 
cept for turbulent entrainment for reasons which become clear 
later. Finally, a summary of our findings is given in Sect. 6. 



2. Physical Conditions and Initial Data 

Our initial helium core models (Tab. 1) with single and double 
convection zones (models M and SC, respectively) are obtained 
from ID stellar evolutionary calculations of a Pop I star (Z=0.02) 
with a mass of 1.25 M Q , and a Pop III star (Z=0) with a mass 
of 0.85 M Q , 4 , respectively. Both models are characterized by 
an off-center helium ignition which results in convection zones 
characterized by a temperature gradient close to the adiabatic 5 
one above the helium burning source. 

The helium cores of both models are composed of a gas 
which is almost completely ionized, as the ionization potentials 
of both He and He + (E B ~ 24.7 eV and E B ~ 54.4 eV, respec- 
tively) are very small compared to the thermal energy, i.e., 



their typical thermal energy is small, 



< 10~ 



(1) 



where kg is Boltzmann's constant, and T > 4.7 x 10 7 K (T > 
1.7 x 10 7 K) holds for the temperature inside the helium core 
of model M (SC) which has a radius r ~ 1.9 x 10 9 cm (r ~ 
5.4 x 10 9 cm). 

In the central part of the models (beneath the convection 
zones) the electron density is so high that the gas is highly de- 
generate. On the other hand, the density of electrons in the sin- 
gle and double convection zone is much lower due to a strong 
expansion that occurred a little earlier in the evolution. Thus, the 
degeneracy has been lifted in the convection zones, i.e., the ratio 
of the Fermi energy Efe of the electrons (Weiss et al. 2004) and 
their typical thermal energy is small, 



Ef.e 



2.47 (0.14) 



(2) 



where E F , E ~ 4.1 x 10~ 8 erg (£ F>£ ~ 1.2 x 10~ 9 erg ) in the 
middle of the convection zone of model M (SC) at a radius r ~ 
7 x10 s cm (r~3.3xl0 9 cm). 

The ions can be described as an ideal, non-relativistic 
Boltzmann gas, because the ratio of their Fermi energy Efj and 



4 Metal-free stars with masses > 1 M Q ignite helium at the center be- 
fore electrons become degenerate (Fujimoto et al. 2000). 

5 The entropy 5 (Fig. 2) and the degeneracy parameter if/ remain al- 
most constant in the convection zones, which is a result of the almost 
adiabatic temperature gradient, i.e., the temperature relation T cc p r ~' 
with the adiabatic exponent y ~ 5/3 holds. Since pT~ 3/2 = f(i/f), the 
degeneracy parameter if/ is constant. 



Efj 
kef 



2 x 10~ 4 (1 x 10~ 5 ) 



(3) 



where E F j ~ 3.4 x 10~ 12 erg {E FJ ~ 1 x 10~ 13 erg ) in the 
middle of the convection zone of model M (SC) at a radius 
r ~ 7 x 10 8 cm (r ~ 3.3 x 10 9 cm). Coulomb forces between 
the ions are negligible, too, as the Coulomb energy correspond- 
ing to the mean ion-ion distance is less than 30% (10%) of the 
thermal energy of the ions for model M (SC) in the middle of the 
convection zone. 

Convection may become turbulent showing random spatial 
and temporal fluctuations. This can be characterized by the di- 
mensionless Reynolds number R e (Landau & Lifshitz 1966) 
which is basically the ratio of inertial to viscous forces. The tur- 
bulent regime is entered once R e exceeds a certain critical value 
Rent, typically being of the order of 10 3 , at which small fluctu- 
ations in the flow are no longer damped. We estimate that the 
Reynolds numbers in the central convection zones of our models 
are 



R e . 



v ■ I ■ p 

n 



10 14 (10 14 ) 



(4) 



where p ~ 10 5 g cm 3 (p ~ 10 3 g cm 3 ), / ~ 10 8 cm 
(/ ~ 10 9 cm), v ~ 10 6 cm s _1 (v ~ 10 5 cm s" 1 ), and 77 ~ 
10 5 g crrT 1 s~' (77 ~ 10 3 g crrr 1 s~') are the typical densities, 
lengths, velocities, and viscosities 6 of the convective flow in 
model M (SC) as predicted by stellar evolutionary calculations. 
These values imply that the flow is highly turbulent, which leads 
to complications when trying to simulate such flows, as turbu- 
lence is an intrinsically three-dimensional phenomenon involv- 
ing a large range of spatial and temporal scales. We recall that, 
in three-dimensional turbulent flow, large structures are unsta- 
ble and cascade into smaller vortices down to molecular scales 
where the kinetic energy of the flow is eventually dissipated into 
heat. 

If L is the largest (integral) scale characterizing a flow, and I 
the scale where viscous dissipation dominates, one has the well 
known relation: 



- ~ R 3/4 



(5) 



In the convection zone, where the Reynolds number R e ~ 10 14 , 
one finds L/l ~ 10 105 . Therefore, the number of grid points /V 
that a numerical simulation would require to resolve all relevant 
length scales is N ~ (L/l) 3 ~ R 9 J A ~ 10 3L5 , which is roughly 22 
orders of magnitude larger than the highest resolutions that can 
be handled by present day computers. 

To account for turbulence on the numerically unresolved 
scales, one usually adopts sub-grid scale models e.g., the quite 
popular one by Smagorinsky (1963), which describe the en- 
ergy transfer from the smallest numerically resolved turbulent 
elements to those at the dissipation length scale using various 
(phenomenological and/or physical) model and flow dependent 
parameters. Here we have not adopted a sub-grid scale model, 
however it has been shown that if one does or does not use such 
a model it seems not to lead to qualitative differences in the hy- 
drodynamic behavior of the core helium flash (Achatz 1995). 



6 The estimate of 77 for a strongly degenerate and completely ionized 
helium gas is based on the formula of Itoh et al. (1987). 
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Table 1. Some properties of the initial models M & SC: Total mass M, stellar population, metal content Z, mass Mhc and radius Rn e 
of the helium core (X( 4 He) > 0.98), ratio of the binding energy of the electrons Eb in the helium core and of their thermal energy 
ksT, ratio of the Coulomb and thermal energy of the ions F cm ,, electron degeneracy parameter in the convection zone W cnv , nuclear 
energy production rate due to helium and hydrogen Lu e burning, temperature maximum T max , and radius r max and density p max 
at the temperature maximum. 
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Fig. 2. Left: Temperature distribution in the helium core in model M (long-dashed), and in model SC (solid) with its stabilized 
counterpart (dash-dotted red), respectively. The two parts of the double convection zone present in model SC are denoted by CVZ-1 
and CVZ-2, respectively. Right: Entropy distribution of model M (solid) and model SC (long-dashed), respectively. 



2. 1 . Initial stellar model M 

The initial model M (Tab. 1, Fig. 2) was obtained with the stel- 
lar evolution code GARSTEC (Weiss & Schlattl 2000, 2007) 
by Achim Weiss, and represents a 1 .25 M Q star at the peak of 
the core helium flash characterized by an off-center temperature 
maximum at the base of a single convection zone sustained by 
helium burning. Additional information about the model can be 
found in Mocak et al. (2008, 2009). 

As we are interested here only in the hydrodynamic evolu- 
tion of the helium core 7 , we consider of model M only the shell 
between r = 2 x 10 8 cm to r = 1.2 x 10 9 cm, which contains the 
single convection zone powered by triple-a burning. Originally, 
the convection zone reaches from 4.72 x 10 8 cm (local pressure 
scale height 2.9 x 10 8 cm) to 9.2 x 10 8 cm (local pressure scale 
height 1.4 x 10 8 cm). From the bottom to the top of the convec- 
tion zone the pressure changes by ~ 1 order of magnitude, from 
6.6 x 10 21 dyn cirT 2 to 7. 1 x 10 20 dyn cirT 2 . We note that both the 
base and the top of the convection zone are located sufficiently 
far away from the (radial) grid boundaries. 

Model M contains the chemical species 'H, 3 He, 4 He, 12 C, 
13 C, 14 N, I5 N, ie O , n O, 24 Mg, and 28 Si. However, since we are 
not interested in details of its chemical evolution, we considered 



7 The helium core is basically a white dwarf sitting inside a red giant 
star. It has a relatively small radius - comparable to that of the Earth - 
but contains almost half of the total mass of the star (Tab. 1). 



only the abundances of 4 He, 12 C, and ie O in our hydrodynamic 
simulations. This is justified as the triple a reaction dominates 
the nuclear energy production rate during the core helium flash. 
For the remaining composition we assume that it can be repre- 
sented by a gas with a mean molecular weight equal to that of 
20 Ne, as its nucleon number agrees well with the average nu- 
cleon number of the neglected nuclear species. 

The stellar model had to be relaxed into hydrostatic equilib- 
rium after it was mapped to the numerical grid of our hydro- 
dynamics code. This was achieved with an iterative procedure, 
which keeps the density distribution of the model almost con- 
stant, but modifies its pressure distribution to achieve hydrostatic 
equilibrium (Mocak 2009). This mapping process has a negligi- 
ble effect on the stellar structure. 

2.2. Initial stellar model SC 

The initial model SC (Tab. 1, Fig. 1 to 3) was computed by 
Simon W. Campbell using the Monash/Mount Stromlo stel- 
lar evolution code (MONSTAR) (Campbell & Lattanzio 2008; 
Wood & Zarro 1981). It corresponds to a metal-free Pop III star 
with a mass of 0.85 M Q near the peak of the core helium flash. 
Metal-free stars with masses > 1 M Q do not undergo the core he- 
lium flash (as opposed to M < 2.2 M Q at solar metallicity). The 
helium core flash commences with a very off-center ignition of 
helium in a relatively dense environment under degenerate con- 
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Fig. 3. Left: Chemical composition of the helium core in model heflpopIH.2d.2 (SC). Right: Nuclear energy production rate as 
a function of radius r. Initial rates (at t=0) are indicated by dotted-black curves. Rates in model heflpopI.3d (SC) at t — 6400 s 
(solid-red), and in model M at t ~ 10 5 s (solid-blue), respectively. 



ditions, and results in a fast growing convection zone powered 
by helium burning that relatively quickly reaches the surround- 
ing hydrogen shell (Fujimoto et al. 1990). This causes sudden 
mixing of protons down into the hot helium core (Fig. 3), and 
leads to rapid nuclear burning via the CNO cycle i.e., a hydro- 
gen flash. Since the core helium flash is still ongoing we refer to 
this as a "Dual Core Flash" (DCF). We note that this event has 
also been referred to as "helium flash induced mixing" (Schlattl 
et al. 2001; Cassisi et al. 2003; Weiss et al. 2004), and "helium 
flash-driven deep mixing" (Suda et al. 2004). The CNO burning 
leads to an increase of the temperature inside the helium-burning 
driven convection zone, and causes it to split into two. The result 
is a lower convection zone still powered by helium burning and 
a second one powered by the CNO cycle (Fig. 3, [right panel]). 
In the following we will refer to the split convection zone as a 
double convection zone. 

The electron degeneracy in the double convection zone has 

already been significantly lifted to tfr 2. It can be shown that 

for a degeneracy parameter tfr < -2, the gas pressure is essen- 
tially that of a nondegenerate gas (Clayton 1968). This confirms 
our previous conclusion based on the ratio of the Fermi and ther- 
mal energy of the electrons (Sect. 2). 

In our hydrodynamic simulations we considered a shell from 
model SC which extends from r — 5 x 10 8 cm to r — 6 x 10 9 cm 
containing the double convection zone. Initially, the inner con- 
vection zone (powered by triple-a burning) covers a region from 
r » 11 X 10 8 cm (local pressure scale height 6 x 10 8 cm) to 
r ~ 35 x 10 8 cm (local pressure scale height 7 x 10 8 cm), while 
the outer convection zone stretches from there up to 54 x 10 8 cm 
(local pressure scale height 5.2 x 10 8 cm). From the bottom 
to the top of the double convection zone the pressure changes 
by ~ 3 orders of magnitude from 1 x 10 20 dyn cirT 2 erg to 
1.4 x 10 17 dyn crrT 2 . Again, we have ensured that the region of 
interest was located sufficiently far away from the radial grid 
boundaries. 

Our hydrodynamic simulations were performed adopting the 
mass fractions of all the species used in the corresponding stellar 
evolutionary calculations, namely 'H, 3 He, 4 He, 12 C, 14 N, and 



16 0. Since the evolutionary model did not include 13 C and 13 N, 
we determined their mass fractions assuming that the CNO cycle 
had been operating in equilibrium. The remaining composition 
is represented by a gas with the molecular weight of 20 Ne. 

The model was relaxed to hydrostatic equilibrium in the 
same manner as in case of model M. This process resulted in 
small fluctuations in the temperature profile (Fig. 2), which were 
smeared out after the onset of convection. 



3. Input physics 

The input physics of our hydrodynamic simulations is identi- 
cal to that one described in Mocak et al. (2008), except for 
the number of nuclear species employed in the simulations 
based on the initial model SC. We use an equation of state that 
includes contributions from radiation, ideal Boltzmann gases, 
and an electron-positron component (Timmes & Swesty 2000). 
Thermal transport was neglected as the maximum amount of 
energy transported by radiation and heat conduction is smaller 
than the convective flux by at least seven (three) orders of mag- 
nitude in model M (SC). Neutrinos act as a nuclear energy 
sink, but carry away less than < 10 2 erg g -1 s _1 . This is a neg- 
ligible amount (especially for the timescales covered by our 
simulations) when compared to the maximum nuclear energy 
production e, which is ~ 10 11 erg g -1 s _1 for model M, and 
~ 10 10 erg g -1 s _1 in model SC, respectively (Fig. 3). 

3.1. Nuclear reactions 

We employed two different nuclear networks for our simula- 
tions, as the nuclear species considered in models M and SC 
differ. The nuclear reaction network used in the hydrodynamic 
simulation based on the initial model M (Tab. 1) consists of four 
nuclei (Sect. 2.1) coupled by seven reactions. The network is 
identical to that one described by Mocak et al. (2008) i.e., 
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The nuclear reactions considered in the hydrodynamic sim- 
ulations based on the initial model SC (Tab. 1) are described by 
a reaction network consisting of nine nuclei (Sect. 2.2) coupled 
by the following 16 reactions: 
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This network reproduces the nuclear energy generation rate of 
the original stellar model very well (Fig. 3). 

Note, that although the value of the temperature maximum, 
T ma x, is higher in model SC than in model M, the energy genera- 
tion rate is lower at T max in model SC, because the helium mass 
fraction X( 4 He) is smaller in that model (0.956 as compared to 
0.970 for model M). 

4. Hydrodynamic code and computational setup 

We use the hydrodynamics code Herakles (Kifonidis et al. 2003, 
2006; Mocak et al. 2008, 2009) which solves the Euler equa- 
tions coupled with source terms corresponding to gravity and nu- 
clear burning. The hydrodynamic equations are integrated with 
the piecewise parabolic method of Colella & Woodward (1984) 
and a Riemann solver for real gases according to Colella & Glaz 
(1984). The evolution of the nuclear species is described by a 
set of additional continuity equations (Plewa & Miiller 1999). 
Self-gravity is implemented according to Miiller & Steimnetz 
(1995) and nuclear burning is treated by the semi-implicit Bader- 
Deufihard scheme (Press et al. 1992). 

We performed one 3D simulation based on the initial model 
M, which covered roughly 27 hrs of stellar evolution (Tab. 2). 
This model (henceforth hefipopI.3d) was evolved on a compu- 
tational grid consisting of a 30°-wide wedge in both 9 and in- 
direction centered at the equator. The small angular size of the 
grid allowed us to achieve a relatively high angular resolution 
(1°) with a modest number of angular zones (N^Ne = 30). 

In addition, we performed two 2D (henceforth models 
heflpopIII.2d. 1 and heflpopHI.2d.2) and one 3D simulation 
(henceforth model heflpopIII.3d) based on the initial model SC 
covering about 1.8 hrs and 0.39 hrs of stellar evolution, respec- 
tively (Tab. 4). We used a computational grid consisting of a 50°- 
wide angular wedge centered at the equator in case of models 
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Fig. 4. Radial velocity distributions for the 3D model 
heflpopI.3d. The dotted and green dashed lines show the time 
(from 10000 s to 30000 s) and angle-averaged radial velocity, 
v r , and velocity modulus, v a b S , respectively. The red dashed line 
shows again the latter velocity, but time-averaged from 80 000 s 
to 99 000 s. The velocity predicted by the mixing-length theory 
(v n ,i,) for the initial model M is shown by the solid blue line. 
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Fig. 5. Convective and kinetic energy fluxes (Fc and Fk, re- 
spectively) as a function of radius averaged (from 33 000 s to 
53 000 s) over about 20 convective turnover timescales for the 
3D model heflpopI.3d. The dotted vertical lines mark the edges 
of the single convection zone in the initial model M according to 
the Schwarzschild criterion. 



heflpopIII.2d. 1 and heflpopIII.3d, and of a 120° wedge in case 
of model heflpopIII.2d.2. 

We imposed reflective boundary conditions in the radial di- 
rection and transmissive ones in the angular directions in all our 
multi-dimensional simulations. 
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Table 2. Some properties of the 3D simulation based on model M: number of grid zones, wedge size w, grid resolution in r, 6 and 
<p direction, estimated Reynolds number R e (Porter & Woodward 1994), characteristic velocity v c of the convective flow, typical 
convective turnover timescale t c „,,, and maximum evolutionary time t„ mx . 
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Fig. 6. Radial distribution of the adiabatic temperature gradient 
V a d (dotted) obtained using the EOS, and of the temperature gra- 
dient of the 3D model heflpopI.3d averaged over the first 460 s 
of its evolution (dashed-red), and over a period of 13000 s be- 
tween 33000 s and 46000 s (solid-black), respectively. The latter 
gradients shown are actually linear fits to the model data. The 
gray shaded region marks the single convection zone CVZ. 



5. Results 

In this section, we first present the characteristics of the hydro- 
dynamic simulation based on the initial model M, i.e., model 
heflpopI.3d, which shows a fast growth of the convection zone 
(Sect. 5.1). A growth rate of this magnitude likely leads to a hy- 
drogen injection phase (Sect. 5.2), which may resemble the one 
seen in the initial model SC, whose hydrodynamic properties are 
discussed in Sect.5.3. 

5.1. Single flash 

Table 2 provides some characteristic parameters of our 3D sim- 
ulation heflpopI.3d based on the initial model M. After convec- 
tion reaches a quasi steady-state in this model, the maximum 
temperature rises at the rate of 80 K s _1 , i.e., only 20% slower 
than predicted by canonical stellar evolution theory. This cor- 
responds to an increase of the nuclear energy production rate 
from 2.4 x 10 42 erg s _1 at t ~ 20000 s to 5.5 x 10 42 erg s _1 at 
t ~ 99000 s. Consequently, the maximum convective velocities 
rise by 26% from about 10 6 cm s _1 to 1.26 x 10 6 cm s _1 . As il- 
lustrated in Fig. 4 these velocities match those predicted by the 
mixing length theory quite well. During the first third of the sim- 
ulation (up to about 30 000 s) the angle and time averaged radial 
velocity in the convective layer exceeds the velocity predicted 
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Fig. 7. Radial distribution of the expansion velocity, v exp , in the 
3D model heflpopl.3d (solid) compared with the expansion ve- 
locity predicted by the stellar evolutionary calculations (dashed- 
blue) for the initial model M. 

for inital model M by the mixing-length theory, v m /,, by about 
20%, while the velocity modulus, v a b s , is about 30% larger than 
v m j r . Towards the end of our simulation the angle and time aver- 
aged modulus of the velocity is about twice as large as v m / f . 

Contrary to our previous study (Mocak et al. 2009), we do 
not find a sub-adiabatic gradient in the outer part of the con- 
vection zone. 8 . The reason for this difference is probably the 
increased grid resolution of our present 3D simulation, which 
results in less heat diffusion due to numerical dissipation, and 
hence a super-adiabatic temperature gradient similar to the ini- 
tial one (Fig. 6) 

At a radius r ~ 5.5 x 10 s cm convection transports al- 
most 90% of the liberated nuclear energy, i.e., 2.1 x 10 42 erg s~' 
(Fig. 5). The energy flux due to thermal transport is negligible 
(Sect. 3). 

We observe internal gravity waves 9 or g-modes in the con- 
vectively stable layers. These g-modes are strongly instigated 
only during certain evolutionary phases (g-mode events) because 
of the intermittent nature of the convective flow. The g-mode 
events are correlated with outbursts in kinetic energy of the con- 

8 A sub-adiabatic gradient does not imply that convective blobs are 
cooler than their environment and that consequently convection ceases 
(the latter is only true when assuming adiabatically rising blobs). It only 
means that blobs cool faster than their environment. 

9 In a convectively stable region, any displaced mass element is 
pushed back by the buoyancy force. On its way back to its original po- 
sition, the blob gains momentum and therefore starts to oscillate. These 
oscillations are called internal gravity waves (Dalsgaard 2003). 
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erg g -1 ) in the 3D model hefipopI.3d. The growth of the size of the convection zone due to the turbulent entrainment, mainly at its 
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vection zone (Meakin & Arnett 2007). It appears that the fre- 
quency of g-mode events increases as the flow gains strength 
(Fig. 9). 

We observe a growth of the size of the single convection zone 
due to turbulent entrainment at the convective boundary on a dy- 
namic timescale (Figs. 8 to 10), and estimate the corresponding 
entrainment speeds by adopting the prescription of Meakin & 
Arnett (2007). Turbulent entrainment involves mass entrainment 
(Fig. 8) rather than a diffusion process, which acts to reduce the 
buoyancy jump at the convective boundary allowing matter to be 
mixed further. The entrainment velocity or the interface migra- 
tion velocity, u e , is given by (see Eq. (32) of Meakin & Arnett 
(2007) 



where Aq is the buoyancy jump, i.e., the variation of the buoy- 
ancy flux q across the convective boundary of thickness h, and 
N 2 the square of the Brunt-Vaisala buoyancy frequency. The 
buoyancy flux is defined as q = gp'v' /p, where g,p, and v are 
the gravitational acceleration, the density, and the velocity, re- 
spectively. The quantities p' and v' denote fluctuations around 



the corresponding mean quantities po and Vo, i.e., p — po + p' 
and v = vq + v', respectively. The thickness h of the convective 
boundary is defined as the half width of the peak in the mass frac- 
tion gradient of 12 C, which varies rapidly at the boundary. The 
entrainment velocities (averaged over 20 convective turnovers at 
the end of the simulation) obtained by this formula are summa- 
rized in Tab. 3 together with those derived from measuring the 
position of the convective boundary in our hydrodynamic simu- 
lations. 

The entrainment velocity derived for our models, v e , are 
calculated by measuring the radial position of the convective 
boundaries defined by the condition X( 12 C) = 2 x 10~ 3 , as 12 C 
is much less abundant outside the convection zone (Fig. 10). The 
expansion velocity (zero in hydrostatic equilibrium) is given by 
v e xp — m r IAnr 2 p, where m r is the partial time derivative of the 
mass m r of a shell of density p at a radius r. 

We find in our models that v e , which includes the expansion 
velocity v exp , agrees very well with the velocity u e predicted by 
theoretical considerations of the entrainment process (see Eq. 6), 
despite the crude estimate of the divergence of the buoyancy flux 
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Fig. 10. Angular averaged 12 C mass fraction as a function of radius near the inner (left) and outer edge (right) of the convection 
zone in the 3D model hefipopI.3d at t — 100000 s. The thick line gives the corresponding temperature stratification, and the vertical 
dotted lines mark the edges of the convection zone at t — s. The observed entrainment velocities v e are given in both panels, too. 



Table 3. Some quantities characterizing the convective bound- 
aries of the 3D model hefipopI.3d : buoyancy jump Aq (in 
erg g -1 s _1 ), square of the Brunt-Vaisala buoyancy frequency 
N 2 , interface migration velocity calculated according to Eq. 6, 
u e , and obtained from our simulation, v e , and expansion velocity 
V exp , at the inner and outer convective interface of thickness h, 
respectively. 
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through Aq/h. The velocity given by the difference v e - v exp , 
which is the actual entrainment speed in our models, differs from 
the theoretically estimate u e by 0.9 m s _1 and 2.0 m s _1 at the 
inner and outer convection boundary, respectively. 

It seems that turbulent entrainment is a robust process which 
has been seen to operate under various conditions in different 
stars (Mocak et al. 2009; Meakin & Arnett 2007). This process 
may be behind the observed Al/Mg anti-correlation (Shetrone 
1996a,b; Yong et al. 2006), which could result from an injection 
of hydrogen into the helium core and a subsequent dredge-up 
(Langer et al. 1993; Langer & Hoffman 1995; Langer et al. 1997; 
Fujimoto et al. 1999). 

5.2. From the single to the dual flash 

If the radial position of the innermost edge of the hydrogen-rich 
layers was fixed at its initial value in model M at r — 1.9X10 9 cm 
(no expansion), the outer convection boundary would reach the 
hydrogen-rich layer due to the turbulent entrainment within only 
17 days, and the helium core would experience a dual core flash 
(DCF) known from low-mass Pop III stars. 

10 d,b = div(q) ~ u e d r b = u e N 2 is a rather crude approximation which 
provides an order of magnitude estimate only (Meakin & Arnett 2007). 



However, the hydrogen layer will initially expand outwards 
at a faster rate than the outer convective boundary (Fig. 7). This 
delays the expected onset of the hydrogen injection a little, as 
the outer convection boundary has to catch up with the hydro- 
gen layer expanding away from the HeCZ. As the expansion ve- 
locities of our hydrodynamic models are biased by the imposed 
reflective boundaries in radial direction, we did not use these val- 
ues here. To get an estimate for the onset of the hydrogen injec- 
tion we instead used the expansion speed of the helium core as 
predicted by stellar evolutionary calculations (Tab. 3), and find 
that the injection of hydrogen into the helium core should take 
place within 23 days. 

This is still a very short time, and even if our estimate for 
the time until the onset of the hydrogen injection was wrong by 
5 orders of magnitude, the injection will occur before the first 
subsequent mini helium flash 1 1 takes place, which is supposed 
to occur in roughly 10 5 years. 

Figure 1 1 shows snapshots from 2D hydrodynamic simula- 
tions (Mocak et al. 2009) of a (single) core helium flash based on 
initial model M (left panel), and of a dual core (helium + hydro- 
gen) flash based on initial model SC (middle and right panels). 
The figure suggests that that even if the core helium flash starts 
out with a single convection zone in a low-mass Pop I star, this 
convection zone may evolve due to its growth by turbulent en- 
trainment (indicated by the arrow in the left panel of the figure) 
into a double convection zone like that found in models of Pop 
III stars. Although the upper boundary of the single convection 
zone is still far away from the hydrogen-rich layers at the end 
of the simulation, it should get there eventually. We find no rea- 



11 When the helium burning shell of the first major helium flash be- 
comes too narrow, it is not able to lift the overlying mass layers. It ex- 
pands slowly, i.e., Sp/p < 0, but remains almost in hydrostatic equilib- 
rium SP/P ~ 0, which in turn leads to a rise of its temperature ST/T > 
(Kippenhahn & Weigert 1990). Hence, helium is re-ignited, but less vi- 
olently than in the first main helium flash. In general, one refers to this 
event as a thermal pulse, but we prefer to call it a mini helium flash. This 
process repeats itself several times until the star settles on the horizontal 
branch. 
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Table 4. Some properties of the 2D and 3D hydrodynamic models based on initial model SC: wedge size w, number of grid points 
in r, 6, and <p direction, estimated Reynolds numbers R e \, R e i (Porter & Woodward 1994), characteristic flow velocities, v c i and v C 2, 
and typical convective turnover timescales, T cm \ and T cnv 2, for the bottom and top part of the convection zone, respectively. t max is 
the maximum evolutionary time. 
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Fig. 11. Snapshots of the spatial distribution of the velocity modulus |v| (in units of 10 6 cm s" 1 ) of a typical 2D hydrodynamic model 
with a single convection zone (left), and for the 2D hydrodynamic model heflpopIII.2d.2 at 1250 s (middle). The right hand panel 
shows the double convection zone highlighted by using the hydrogen contrast AX('H) = 100 x (X^H) - <X( 1 H)) e )/<X( 1 H)) e at the 
same time, where ()g denotes a horizontal average at a given radius. The arrow indicates the growth of the single convection zone, 
while the dotted line represents the border between the helium and hydrogen rich layers. 



son in any of our 2D or 3D simulations including the new model 
heflpopI.3d why turbulent entrainment should cease before the 
outer convective boundary reaches the edge of the helium core. 
Actually, as the maximum temperature of the helium core grows 
and the convective flow becomes faster, entrainment will even- 
tually speed up (Mocak et al. 2009). 

On the other hand, subsequent mini helium flashes are un- 
likely to occur, because we estimate that at the observed entrain- 
ment velocity the inner convective boundary will reach the cen- 
ter of the star in > 90 days. Note in this respect that the fast en- 
trainment speed of the inner convective boundary derived from 
our previous 2D models (Mocak et al. 2008) is not confirmed by 
the new 3D model. The entrainment rate is actually slower in the 
new 3D model, which was expected (Mocak et al. 2009). 

The turbulent entrainment at the inner convection boundary 
heats the layers beneath at a rate 6T/6t ~ 630 K s _1 , i.e., it is 
quite efficient in lifting the electron degenaracy in the matter be- 
low the convection zone. 

According to the above discussion we propose a somewhat 
speculative scenario for the core helium flash in low-mass Pop I 
stars. These stars ignite helium burning under degenerate condi- 



tions and develop a single convection zone, which at some point 
extends due to turbulent enrainment up to their hydrogen-rich 
surface layers. The convection zone eventually penetrates these 
layers, and dredges down hydrogen into the helium core. This ig- 
nites a secondary flash driven by CNO burning, which together 
with triple-a burning and inwards turbulent entrainment leads 
eventually to the lifting of the core's degeneracy, i.e., the star 
will settle down on the horizontal branch. 

5.3. Dual flash 

We have performed three simulations of the core helium flash 
based on the Pop III initial model SC which possesses two con- 
vection zones sustained by helium and CNO burning, respec- 
tively. Some characteristic properties of these dual flash simula- 
tions are summarized in Table 4. 

5.3.1. Model heflpoplll.2d.2 

Despite a nuclear energy production rate due to CNO burning in 
the outer convection zone (e max ~ 1.4 x 10 10 erg g" 1 s _1 , which 
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Fig. 12. Temporal evolution of the radial distribution of the 
(color coded) logarithm of the angular averaged radial (v^/2; 
upper panel) and angular Vg/2; lower panel) component of 
the specific kinetic energy (in erg g~') of the 2D model 
heflpopIII.2d.2. 



is roughly eight times higher than that due to triple-a burning in 
the inner zone e max ~ 1.7 x 10 9 erg g -1 s -1 ), convective motions 
first appear within the inner convection zone after about 200 s. 
The onset of convection in the outer zone is delayed until about 
500 s (Fig. 12, Fig. 16). After some time the model relaxes into a 
quasi steady-state, where the r.m.s values of the angular velocity 
are comparable or even larger than those of the radial component 
(Fig. 12). This property of the velocity field implies the presence 
of g-modes or internal gravity waves (Asida & Arnett 2000), 
which is a surprising fact. G-modes should not exist in the con- 
vection zone according to the canonical theory, as any density 
perturbation in a convectively unstable zone will depart its place 
of origin exponentially fast (Kippenhahn & Weigert 1990) until 
the flow becomes convectively stable. 

When convection begins to operate in both zones, the total 
energy production rate temporarily drops by ~ 20%, but con- 
tinues to rise at a rate of 2.8 x 10 36 erg s~ 2 after t as 2000 s. 
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Fig. 13. Adiabatic temperature gradient V ac / (solid-black), tem- 
perature gradient of the 2D model heflpopIII.2d.2 averaged from 
1000 s to 3200 s (dashed-red), and temperature gradient of the 
initial model SC (solid-blue), respectively, as a function of ra- 
dius. 



Nevertheless, the convective flow decays fast in both convec- 
tion zones - at a rate of 4 x 10 40 erg s~'(Fig. 16), probably 
because the initial conditions represented by the stabilized ini- 
tial model are too different from those of the original stellar 
model. They disfavor convection since the stabilized model has a 
slightly lower temperature gradient than the original one (Fig. 2, 
Fig. 13). However, we note that stabilization of the initial model 
was essential to prevent it from strongly deviating from hydro- 
static equilibrium. 

Shortly after convection is triggered in both the outer and the 
inner zone this double convection structure vanishes, and after 
about 2000 s there is no evidence left for two separate convection 
zones (Fig. 12). 

Due to the relatively short temporal coverage of the evolu- 
tion and due to the decaying convective flow, we did not ana- 
lyze the energy fluxes and turbulent entrainment of this model. 
Penetrating plumes do not exist in the convection zone (initially 
determined by the Schwarzschild criterion), as it is dominated 
rather by g-modes. Hence, firm conclusions are difficult to de- 
rive, but we plan to address this issue elsewhere. We are now 
going to introduce some basic characteristic of internal gravity 
waves or g-modes that are required to draw further conclusions. 

G-modes In a convectively stable region, any displaced mass 
element or blob (density perturbation) is pushed back by the 
buoyancy force. On its way back to its original position, the blob 
gains momentum and therefore starts to oscillate around its orig- 
inal position. Assuming, the element is displaced by a distance 
Ar, has an excess density Ap, and is in pressure equilibrium with 
the surrounding gas (AP = 0), one can derive an equation for the 
acceleration of the element: 



d 2 (Ar) 



V ele " V 



surr + — v p 



Ar 



(7) 
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Fig. 14. Radial distributions of the (square of the) Brunt-Vaisala 
buoyancy frequency A' 2 in the 2D model heflpopIII.2d.2 aver- 
aged between 1480 s and 6000 s (top), and in the 3D model 
heflpopl.3d averaged between 6600s and 100000 (bottom), re- 
spectively. The angular and temporal variation of A^ 2 at a given 
radius are indicated by the gray shaded region. The inserts show 
a zoom of the region around A^ 2 = to enlarge the variations of 
Af 2 in the convection zone which are.<10~ 3 . 



where g is the gravitational acceleration, H p = \p/d r p\ is the 
pressure scale height, V = iiln T/dlnP, = dln/j/dlnP, 6 = 
31np/<91n T, and <p = <91np/<91n/i. 

Let us assume now that the element, after an initial displace- 
ment Aro, moves adiabatically (V e j e = V a( j) through a convec- 
tively stable layer. The element is accelerated back towards its 
equilibrium position and starts to oscillate around this position 
according to the solution of Eq. (7): 

Ar = Ar e* Nr , (8) 
where the (square of the) Brunt-Vaisala frequency is given by 



^ = |b( V ad-Vsurr + f^) 



(9) 



In a convectively unstable region (assuming = 0), Eq.(9) im- 
plies Af 2 < 0, i.e.,N is imaginary. Thus, the displaced element 
moves exponentially away from its initial position, instead of 
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r (10 s cm) 

Fig. 15. Hydrogen mass fraction as a function of radius for 
the 2D model heflpopIII.2d.2 at t — 6405 s. The vertical lines 
mark the initial border between hydrogen and helium rich lay- 
ers (dashed-red), and the layer (dotted) where the timescale for 
proton capture on 12 C equals 10 2 s (T ~ 1 x 10 s K). 



oscillating around it, as it is the case for g-modes or internal 
gravity waves and A^ 2 > 0. 

G-modes appear in layers of gas stratified under gravity and 
are spatial oscillatory displacements of density perturbations. 
The dispersion relation of such density displacements, assumed 
to vary as exp[/(k • r - o)f)], reads (Dalsgaard 2003) 



N 1 



l+k 2 r /k 2 h 



(10) 



where u is the temporal frequency of the density displacements, 
and k r and k/, are the radial and horizontal components of the 
wave number k, respectively. The dispersion relation tells us 
that any density perturbation must (under influence of buoyancy 
forces) displace matter horizontally to propagate vertically. This 
will give rise to matter motion resembling horizontal fingers as a> 
approaches the Brunt-Vaisala frequency for large k/, (k 2 /k 2 t — > 
and cu 2 — > N 2 ). 

We find such horizontal structures in our models (Fig. 11) 
visible mainly in the radiative layer of the splitted convection 
zone (Fig. 1 1). By decomposition of the specific kinetic energy 
density of the model into the radial (v 2 /2) and horizontal (v?/2) 
component, we also find that the horizontal displacements are 
characterized by higher values of the kinetic energy density 
compared to the corresponding values of the vertical displace- 
ments already 2000 s after the start of the simulation (Fig. 12). 
Additionally, A^ 2 is mostly positive in these models (Fig. 14), in- 
dicating convective stability throughout the double convection 
zone. This proves the existence of internal gravity waves in the 
decaying double convection zone. The situation is different in 
our 3D model heflpopl.3d, where A^ 2 is (on average) small and 
negative everywhere in the single convection zone (Fig. 14). 

Within the double convection zone originally determined by 
the Schwarzschild criterion the temperature gradient drops ev- 
erywhere below the adiabatic one (Fig. 13). It does not imply 
that convection must cease. Even if the temperature gradient is 
not everywhere super-adiabatic nuclear burning may create hot 
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Fig. 16. Left panel: Temporal evolution of the nuclear energy production rate E (left) in units of the solar luminosity L and the 
kinetic energy Ek (right) for the models hefipopIII.3d (diamonds-black), hefipopIII.2d.l (solid-red), and hefipopIII.2d.2 (dashed- 
blue), respectively. The green stars give the behavior of a 3D model with the same properties as model hefipopIII.3d, but with 
nuclear burning switched off. The vertical arrows, labeled 1. and 2., mark the onset of convection in the lower and upper convection 
zone of the double convection structure, respectively. 




Fig. 17. Snapshots of the spatial distribution of the velocity modulus |v| (in units of 10 6 cm s ) for a typical 3D hydrodynamic 
model with a single convection zone (left), and for model hefipopIII.3d at t — 1300 s with its double convection zone (middle) 
in the meridional plane = 0°. The right panel shows the hydrogen convection zone highlighted by using the hydrogen contrast 
AX('H) = 100 x (X( l H) - (X( l H))g) I (X( l H))g, where ()e denotes a horizontal average at a given radius. The arrow indicates the 
growth of the single convection zone, while the dotted line represents the border between the helium and hydrogen rich layers . 



blobs, which although cooling faster than the environment can 
still be hotter than the latter, and thus can rise upwards. 

This supports our conclusions based on the distribution of 
the Brunt-Vaisala frequency, which might be, however, a result 
of insufficient resolution, as the gradient increases with increas- 
ing resolution. Consequently, we do not find the typical 2D con- 



vective pattern characterized by vortices in the double convec- 
tion zone at later times t > 2000 s. 

Radiative barrier Stellar evolutionary calculations of low-mass 
Pop III stars predict that the helium flash-driven convection zone 
splits into two, when the hydrogen-burning luminosity (driven 
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Fig. 18. Temporal evolution of the radial distribution of the (color coded) logarithm of the angular averaged kinetic energy density 
(in erg g' 1 ) of models hefipopHI.2d. 1 (left) and hefipopIII.3d (right), respectively. The horizontal dotted lines mark the boundaries 
of the double convection zone one part being sustained by the helium burning (CVZ-1) and the other one by the CNO cycle (CVZ-2). 



by the entrainment of H) exceeds the helium-burning luminos- 
ity. At this point, a radiative barrier is created between both con- 
vective zones, as energy flows inwards from the layers where 
hydrogen-burning takes place. The radiative barrier is thought 
to prevent the flow of isotopes into the helium burning layers 
and vice versa, hence preventing the reaction 13 C (a, n) 16 to 
become a source of neutrons. Also an eventual mixing of iso- 
topes from the helium-burning layer into the stellar atmosphere 
should be inhibited. However, this scenario is difficult to prove 
due to numerical problems arising when modeling this event in 
ID (Hollowell et al. 1990). 

Our hydrodynamic model heflpopIII.2d.2 with the double 
convective zone shows that the radiative barrier allows for some 
interaction between both zones via g-modes (Fig. 11, Fig. 12). 
In addition, there is some mixing of hydrogen into the radia- 
tive layer, which was initially completely devoid of hydrogen 
(i.e., X('H) = 10 30 ). Hydrogen must have been mixed there 
either by convective motion from the hydrogen-rich layers or 
dredged down by penetrating convective plumes from the lower 
convection zone. The downward mixing of hydrogen extends to 
a radius of ~ 2.84 x 10 9 cm (X('H) ~ 10~ 29 ) by the end of 
our simulation heflpopIII.2d.2 (Fig. 15). It is likely that deeper 
mixing of hydrogen into the helium burning layers is not pos- 
sible, since protons are captured via the reaction 12 C (p, y) 13 N 
on timescales shorter than that on which protons are mixed in- 
wards (Hollowell et al. 1990). At a temperature T ~ 10 8 K 
(r ~ 2.21 x 10 9 cm) the proton lifetime against capture by 12 C is 
as short as ~ 10 2 s (Caughlan & Fowler 1988). This is an order of 
magnitude smaller than the observed initial convective turnover 
timescales (Tab. 4). 

5.3.2. Simulation heflpoplll.2d.1 and heflpoplll.3d 

We now discuss the qualitative behavior of 2D and 3D Pop III 
models which were simulated using the same number of ra- 
dial and angular zones, but which have a lower radial and an- 
gular resolution than the 2D model heflpopIII.2d.2 discussed in 
the previous subsection. The quantitative properties of the con- 



vection zone of these models will obviously be different as an 
increased grid resolution implies less numerical viscosity and 
larger Reynolds numbers. We again stress here that the charac- 
teristic Reynolds numbers in our 2D and 3D simulations (Tab. 4) 
are still many orders of magnitude smaller than the values pre- 
dicted by theory (see Sect. 2). 

The comparison between the 2D and 3D simulations, 
heflpopIII.2d. 1 and heflpopIII.3d, provides important informa- 
tion on the impact of the symmetry restriction imposed in 
the 2D models. Contrary to the 2D models, our 3D hydrody- 
namic simulations of turbulent flow performed with the PPM 
scheme (Sect. 4) are geometrically unconstrained, i.e., in the in- 
ertial regime turbulent eddies can decay along the Kolmogorov 
cascade down to the finest resolved scales (Sitine et al. 2000). 

Due to the large computational cost we evolved the 
3D model heflpopIII.3d and the corresponding 2D model 
heflpopIII.2d.l for 0.39 hrs of stellar life, only. We find the fol- 
lowing qualitative differences between the 3D and 2D model 
(Fig. 16): (i) in 3D convection starts earlier in the outer part of 
the double convection zone, (ii) convective velocities are smaller 
there, and (iii) the convective structures have a plume-like shape 
in 3D (Fig. 17) and are vortex-like in 2D (Fig. 1 1). On the other 
hand, the models also exhibit the following common qualitative 
evolutionary properties (Fig. 16): (i) a decrease of the total nu- 
clear energy production rate, (ii) a decrease of the maximum 
temperature, (iii) a decay of the velocity field in the convection 
zones (Fig. 18), and (iv) the presence of internal gravity waves 
in the radiative barrier. 

The differences observed between the 2D and 3D simulation 
do not come at a surprise, as it is well known that 2D simula- 
tions lead to an overestimate of the flow velocities (Muthsam 
et al. 1995; Meakin & Arnett 2006). On the other hand, the com- 
mon properties of the 2D and 3D model are also shared by the 
high resolution simulation heflpopIII.2d.2, except for the nuclear 
energy production rate, which does not decrease after convec- 
tion is fully established. This implies that both our 2D model 
heflpopIII.2d. 1 and 3D model heflpopIII.3d are not sufficiently 
well resolved, although they show the most important charac- 
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teristics of the high resolution model heflpopIII.2d.2 described 
in Sect. 5.3.1, i.e., the presence of a decaying convective flow 
in both convection zones which are later dominated by inter- 
nal gravity waves. This also holds for the intermediate radiative 
layer. 

Contrary to the low resolution 2D model heflpopIII.2d.l, the 
convective velocities found for the 3D model heflpopIII.3d 12 in 
the inner convection zone sustained by helium burning match 
those predicted by stellar evolutionary calculations relatively 
well, although the modulus of the velocity is about a factor of 
two larger. In the outer part of the convection zone, sustained by 
the CNO cycle, the modulus of the velocity and the individual 
velocity components are smaller by more than a factor of two. 

Interestingly, we find convection to be triggered sponta- 
neously in these simulations - even without nuclear burning. 
This is highlighted by the fact that the temporal evolution of 
the kinetic energy of the 3D model heflpopIII.3d with no nu- 
clear nuclear energy production is almost identical to that of the 
corresponding model with burning switched on (Fig. 16). Thus, 
we conclude that the hydrodynamic convective flow observed in 
our models is mainly driven by the adopted temperature gradi- 
ent which is inherited from the ID stellar model, and is only 
partially sustained by nuclear burning within the hydrodynamic 
simulation. 

6. Summary 

We have performed and analyzed a 3D hydrodynamic simulation 
of a core helium flash near its peak in a Pop I star possessing a 
single convection zone (single flash) sustained by helium burn- 
ing. The simulation covers 27hrs of stellar life, or roughly 100 
convective turnover timescales. In addition, we performed and 
analyzed 2D and 3D simulations of the core helium flash near 
its peak in a Pop III star which has a double convection zone 
(dual flash) sustained by helium and CNO burning, respectively. 
These simulations cover only 1.8 hrs and 0.39 hrs of stellar life, 
respectively, as convection dies out shortly after it appears. 

The convective velocities in our hydrodynamic simulation of 
the single flash model and those predicted by stellar evolutionary 
calculations agree approximately. Contrary to our previous find- 
ings, the temperature gradient in the convection zone remains 
superadiabatic, probably because of the increased spatial resolu- 
tion of these simulations as compared to our old models. As ex- 
pected, the simulation shows that the convection zone grows on 
a dynamic timescale due to turbulent entrainment. This growth 
can lead to hydrogen injection into the helium core as predicted 
by stellar evolutionary calculation of extremely metal-poor or 
metal-free Pop III stars. Hydrogen injection leads to a split of the 
single convection zone into two parts separated by a supposedly 
impenetrable radiative zone. Our hydrodynamic simulations of 
the double convection zone show that the two zones vanish as 
their convective motion decays very fast. However, this result 
may be caused by an insufficient spatial grid resolution or prob- 
ably because the conditions represented by the stabilized initial 
model are a bit different from those of the original stellar model. 
While the convective velocities in our 2D hydrodynamic models 
do not match those predicted by stellar evolutionary calculations 
for the double convection zone at all, a rough agreement is found 
in our 3D model for the velocities in the inner convection zone 
sustained by helium burning. 

12 The convective velocities are measured just after convection ap- 
pears for the first time during the simulation, as the convective flow 
decays very fast later. 
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